Transitions in the horizontal transport of vertically vibrated granular layers 
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Motivated by recent advances in the investigation of fluctuation-driven ratchets and flows in ex- 
cited granular media, we have carried out experimental and simulational studies to explore the 
horizontal transport of granular particles in a vertically vibrated system whose base has a sawtooth- 
shaped profile. The resulting material flow exhibits novel collective behavior, both as a function of 
the number of layers of particles and the driving frequency; in particular, under certain conditions, 
increasing the layer thickness leads to a reversal of the current, while the onset of transport as a 
function of frequency occurs gradually in a manner reminiscent of a phase transition. Our exper- 
imental findings are interpreted here with the help of extensive, event driven Molecular Dynamics 
simulations. In addition to reproducing the experimental results, the simulations revealed that the 
current may be reversed as a function of the driving frequency as well. We also give details about the 
simulations so that similar numerical studies can be carried out in a more straightforward manner 
in the future. 
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I. INTRODUCTION 

The best known and most common transport mech- 
anisms involve gradients of external fields or chemical 
potentials that extend over the distance traveled by the 
moving objects. However, recent theoretical studies have 
shown that there are processes in far from equilibrium 
systems possessing vectorial symmetry that can bias 
thermal noise type fluctuations and induce macroscopic 
motion on the basis of purely local effects. This mech- 
anism is expected to be essential for the operation of 
molecular combustion motors responsible for many kinds 
of biological motion; it has also been demonstrated ex- 
perimentally in simple physical systems indicating 
that it could lead to new technological developments such 
as nanoscale devices or novel types of particle separa- 
tors. Motivated by both of these possibilities, as well as 
by interesting new results for flows in excited granular 
materials |3j-|8|], we have carried out a series of experi- 
mental and simulational studies that explore the manner 
in which granular particles are horizontally transported 
by means of vertical vibration. 

In the corresponding theoretical models - known as 
"thermal ratchets" - fluctuation-driven transport phe- 
nomena can be interpreted in terms of overdamped Brow- 
nian particles moving through a periodic but asymmetric, 
one-dimensional potential in the presence of nonequilib- 
rium fluctuations |p[-p^|. Typically, a sawtooth-shaped 
potential is considered, and the nonlinear fluctuations 
are represented either by additional random forces or by 
switching between two different potentials. Collective ef- 
fects occurring during the fluctuation-driven motion have 
also been considered |l3|-|l5| , leading to a number of un- 
usual effects that include current reversal as a function of 
particle density. Here we investigate an analogous trans- 
port mechanism for granular materials. By carrying out 
experiments - both real and numerical - on granular ma- 



terials vibrated vertically by a base with a sawtooth pro- 
file, it is possible to achieve a fascinating combination 
of two topics of considerable current interest - ratchets 
and granular flows. A number of recent papers have fo- 
cused on vibration-driven granular flow, and the details 
of the resulting convection patterns have been examined, 
both by direct observation j|,|]|7]|| and by magnetic reso- 
nance imaging |fjjf6f| . Granular convection has also been 
simulated numerically by several groups; the study most 
closely related to the present work deals with the hori- 
zontal transport that occurs when the base is forced to 
vibrate in an asymmetric manner fllTf. 



II. EXPERIMENTS 

We have experimentally investigated of the horizontal 
flow of granular material confined between two upright 
concentric cylinders undergoing vertical vibration. In or- 
der to induce transport, the height of the annular base be- 
tween the cylinders has a periodic, piecewise-linear profile 
(in other words, it is sawtooth- like) . The experimental 
setup is also described in Ref. |f8) , and is briefly reviewed 
here. 



A. Apparatus 

Figure |l| shows a schematic view of the experimental 
apparatus. To achieve a quasi- two-dimensional system 
without boundaries in the direction of the expected flow 
the granular material is placed between two concentric 
glass cylinders . The mean diameter of the cylinders is 
9.2 cm, while the gap between the cylinders is 3.5 mm. A 
ring filling the gap between the cylinders, with a sawtooth 
profile on its upper surface, is mounted on the base of the 
container; the ring is made of either PVC (soft plastic), 
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hard plastic or aluminum, and different sawtooth shapes 
are used. The entire assembly is vertically vibrated with 
a displacement that depends sinusoidally on time. 



B. Results 




FIG. 1. The schematic view of the experimental apparatus. 




FIG. 2. The average horizontal velocity of the beads as a 
function of height of the granular layer. Here we characterize 
the height of the layer by n = N/k, where N is the total num- 
ber of particles and k is the number of particles in a single 
layer. The various curves represent measurements for vari- 
ous sawtooth shapes (given by their width w, height h, and 
asymmetry parameter a, which is the ratio of the horizontal 
projection of the left part of a tooth to its total extension) 
and materials. The shape of the sawtooth can be seen in the 
figure, the filled sawtooth represents hard plastic material, 
the others are made of P VC (soft plastic) . The particles are 
glass balls. The amplitude and the frequency are A — 2 mm, 
/ = 25 Hz. 

Monodisperse glass balls are used in the experiments, 
which are nearly spherical with diameter 3.3 mm ±2%. 
As shown in the inset of Fig. g, the size of the sawtooth 
is in the same range as that of the particles. 



Provided the frequency is sufficiently large, the verti- 
cal vibration causes horizontal flow of the entire granular 
layer. This bulk motion is reproducible over repeated 
experiments. The average flow velocity is determined by 
tracking individual tracer particles visible through the 
transparent cylinder walls. In order to average out fluctu- 
ations, the particles are allowed to travel large distances; 
depending on the size of the fluctuations this distance is 
between 1.5 and 6 m (equal to 5-20 times the circumfer- 
ence of the system). Each point shown in the graphs is 
an average over 3-6 tracer particles. 

Figure 2 shows the horizontal flow velocity as a func- 
tion of the number of particles for various possible sys- 
tems. The actual sawtooth and particle shapes are also 
indicated. Positive velocities arc defined as follows: Mov- 
ing in the positive direction the first edge of a saw- 
tooth is the steeper one; i.e., from left to the right for 
these cases. The vibration amplitude and frequency are 
A = 2 mm and / = 25 Hz; the dimensionless acceleration 
r = (2irf) 2 A/g is an important quantity for vibrated 
granular systems, so that here T = 5. 

We have observed a variety of qualitatively different 
kinds of behavior, some of which can be interpreted by 
simple geometrical arguments. The most surprising phe- 
nomenon is that in certain cases the velocity changes sign; 
in other words, the flow direction depends on the layer 
thickness. In some cases the curves are monotonically 
decreasing, while others have well defined maxima. Al- 
tering the particle shape reduces the velocity and shifts 
the location of the maximum, but the shape of the curve 
remains unchanged (results not shown here). Likewise, 
changing the elasticity of the base does not alter the 
curves qualitatively. The only feature common to the dif- 
ferent curves is that beyond a certain layer thickness the 
velocity magnitude decreases as further layers are added. 

Figure || shows the T dependence of the flow velocity 
for a system of 200 balls (amounting to 4 layers) for con- 
stant A. Flow occurs only above a critical acceleration 
r c ~ 1.7. Above this critical value the velocity appears 
to follow a power law 



v(T) oc (r - r c ) 



0.48 



(1) 



suggesting that the onset of flow resembles the kind of 
phase transition observed in hydrodynamic instabilities 
such as thermal convection fll9fl . In this particular case 
the exponent we obtained is close to 1/2 indicating that 
the transition is more like a bifurcation. However, for 
other parameters (see, e.g., the simulational results pre- 
sented in section [V ) we obtained other exponents as well. 
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FIG. 3. Horizontal velocity (v x ) as a function of the dimen- 
sionless acceleration V at constant amplitude (A = 2 mm). 
The experiment is for a strongly asymmetric aluminum saw- 
tooth (w = 12 mm, h = 7 mm, a = 0) and 200 glass balls. In 
the inset we display the data on a log-log scale for (v x ) close 
to the transition as a function of T* = (r — r c )/r c where 
T c = 1.7. The slope of the fitted line is 0.48. 



III. EVENT DRIVEN SIMULATIONS 



Simulations have been useful in the studies of granu- 
lar systems pO| , pi| . Here we perform event driven sim- 
ulations of inelastic particles with an additional shear 
friction in a two-dimensional system whose base has a 
sawtooth-shaped profile. In event driven simulations the 
difficult part of the algorithm is determining the next 
event (e.g., the next collision), and the motion of the 
part icles is calculated analytically between two events 
p4]-p6|. The event driven simulations discussed below 
represent an alternative to the more common Molecu- 
lar Dynamics calculations assuming differentiable inter- 
action potentials between the particles p7| . 

The particles are modeled as disks with a given radius 
r,;. These disks can rotate, their moment of inertia is 

— rriirj (m; is the mass of the ith particle). The motion 

of the beads is simulated in two-dimensional cell (see Fig. 
|J). To realize the topology of the experiment, periodical 
boundary condition is applied in the horizontal direction. 
The height of the cell is chosen in such a way that the 
number of particles hitting the upper wall is negligible. 
The base is oscillating sinusoidally with frequency / and 
amplitude A. The geometry of the base can be seen in 
Fig. H in detail. The shape of a tooth can be described 
by three parameters: its height h, width w and asymme- 
try parameter a. This asymmetry parameter is a = l/w, 
where I is the distance of the projection of the sawtooth's 
top point from the left side point of the unit. The top 
point of a tooth is rounded by an arc which smoothly 
joins the two sides. 
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FIG. 4. Visualization of the arrangement used in the sim- 
ulations. The upper wall is fixed, the sawtooth shaped base 
is oscillating sinusoidally. Periodic boundary condition is ap- 
plied in horizontal direction. 




FIG. 5. The geometry of the sawtooth base. The shape of 
a tooth can be described by three parameters: its height h, 
width w, and asymmetry parameter a — l/w, where I is the 
horizontal projection of the left part of the tooth. The top 
of the teeth are rounded by an arc with radius s smoothly 
joining to both sides. 

We consider collisions during which the particles can 
stick to or slide along each other's surface. For calcu- 
lating the velocities of two particles after a collision, a 
simple model is used which corresponds to recent colli- 
sion models f28|-|30| with /3o = 0, i.e., without tangential 
restitution. 

The algorithm described in |26| has been implemented 
in our simulation with some extensions, e.g. gravitational 
field, rotating particles, sinusoidally moving objects. In 
the Appendix we give some details of the calculations 
specific to this simulation so that similar studies could 
be carried out or reproduced in the future in a more di- 
rect manner. 

A typical problem arising in event driven granular sim- 
ulations is called inelastic collapse jnj. One possibility 
to avoid the numerical difficulties associated with it is 
to make the coefficient of normal restitution (e) velocity 
dependent [ pl"| , p2[ . We use a very simple rule in our sim- 
ulation: e is constant until v' n > v™ m (v' n is the relative 
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normal velocity of the colliding particles after the colli- 
sion), however, if v' n would be less than v™ n , it is set to 
this minimal normal velocity: v' n = v™ m . Therefore, the 
relative normal velocity cannot be less than i>™ m between 
two objects after their collision. Using w™ m = 5 mm/s 
is enough to avoid inelastic collapse, and still this extra 
rule does not influence the dynamics of the system, as far 
as the transport is concerned. We checked this by deter- 
mining the proportion of collisions when this rule had to 
be used and checking the change in the average velocity 
as the threshold value was lowered. We found that in 
most situations the number of collisions when this rule 
has to be applied is negligible, and even in dense systems 
we did not detect a change (beyond error bars) in the 
transport velocity when t>™ m was decreased. The simu- 
lation program was written in C programming language, 
and it ran on personal computers with Linux operating 
system. 



increased, the direction of the transport becomes posi- 
tive. We chose n = 1, and the frequency was / = 25 Hz 
(r = 5.0). The width of the cell was 24w = 144 mm, 
therefore one layer contained about 70 particles. These 
simulations lasted for 50 internal seconds, the following 
results represent averages of two runs for each parameter 
sets (with different initial conditions). 

In Fig. H (v x ) can be seen vs. a (the height was fixed 
at four different values). The first thing worth to note is 
that (v x ) tends to zero as the shape becomes symmetrical 
(i. e. a tends to 0.5). The shape of the plots depends on 
the height h of the sawtooth. At h = 3 mm the velocity is 
positive and decreases monotonically, while at h = 5 mm 
it is negative and has a minimum at a ~ 2.5, and when 
the height is 9 mm, it is still negative but monotonically 
increasing. These curves are helpful in understanding the 
mechanism of negative transport. 



IV. RESULTS 

The system we studied is a complex one, with many 
parameters. Since it is impossible to explore the depen- 
dence of its behavior on each of the parameters, some of 
them were fixed in all of the simulations. (Still there re- 
mained quite a few parameters to vary, the results shown 
in this section are selected from simulations with more 
than 5000 different parameter sets.) We fixed two pa- 
rameters: the width of a sawtooth was w — 6 mm and 
the amplitude of the oscillation was A = 2 mm. Accord- 
ing to previous results |33| , the direction of the horizontal 
transport does not depend very much on the friction co- 
efficient /i. It may, however, depend on the coefficient of 
restitution e, but due to the great number of other pa- 
rameters, e was kept at fixed value of 0.8 in all parameter 
sets except for one series when it was varied from 0.4 to 
1.0. For the sake of simplicity, both e and \i was chosen 
to be the same in particle-particle and particle-wall col- 
lisions, although the simulation allows setting different 
values for different type of collisions. The radii of the 
particles varied from 1.05 mm to 1.155 mm uniformly, to 
avoid the formation of a hexagonal structure which often 
appears in two dimensions when the system freezes. The 
masses were equal, hence it did not matter what their 
actual value was. 



A. Dependence on the sawtooth shape 

First we investigated how (v x ) depends on the shape of 
the sawtooth (other parameters were kept fixed), where 
the average has been taken over both time and the in- 
dividual particles. We intended to find sawtooth shapes 
resulting in negative transport. According to the exper- 
iments, negative transport occurs only for a few particle 
layers (n — 0...3), and as the number of particles is 



In Fig. the dependence of the horizontal transport 
on the height of the sawtooth can be seen. The shapes of 
the graphs are similar in case of a — 0.05 and a — 0.333: 
(v x ) is positive if h is small, and decreases and turns 
to negative with increasing h. Finally, it becomes zero 
for high h values. More detailed investigations showed 
that it is true for nearly all possible values of a (except 
for a ~ 0.5, since then the transport vanishes). In con- 
clusion, the horizontal transport can be negative if the 
height of the sawtooth is between 4 and 16 mm and if 
the asymmetry parameter a is less then about 0.4. In 
other cases it is positive or zero (if the sawtooth is sym- 
metric). It is important to note that these results are 
valid only for a small number of particles. As it will be 
shown, increasing the number of particles the negative 
transport vanishes and turns into positive in all cases. 
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FIG. 6. The horizontal transport as a function of the asym- 
metry of the sawtooth. The shape of the curves depends on 
the height h of the sawtooth. 
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FIG. 7. Horizontal transport as a funcion of height of the 
sawtooth. The shape does not depend very much on the asym- 
metry of the sawtooth. 
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FIG. 8. The horizontal transport as a function of layer 
width at three different frequencies, when the shape of the 
sawtooth is: h — 6 mm and a — 0.25. The system width was 
lOw = 60 mm, one run lasted for 30 internal seconds. The 
data is average of 7 runs. 



B. Dependence on the layer width 



The results in the previous section show that if there 
are only few particles, then depending on the shape of 
a sawtooth, the transport can be negative or positive. 
Now we show what happens if more and more particles 
are added into the system. 

In Fig. U one can see that negative transport becomes 
zero and turns to positive as the layer width is increased. 



(Unfortunately the event-driven simulation becomes un- 
feasible as the density increases with increasing particle 
number, therefore n ~ 5 is the maxium layer width we 
can have in our simulation.) If the sawtooth shape is 
such that even for a few particles the transport is posi- 
tive, then with increasing particle number the transport 
decreases (see Fig. ||). 
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FIG. 9. The horizontal transport as a function of layer 
width at three different frequencies. The shape of the saw- 
tooth: h = 3 mm and a = 0.1. The data shown here is average 
of 7 runs, one run lasted for 30 internal seconds. 
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FIG. 10. Horizontal transport as a function of driving fre- 
quency, with six different layer widths, where the relation 
between the frequency and the dimensionless acceleration is 
/ = (rff) 1/2 A- 1/2 (27r)~ 1 , with A = 2 mm. The shape of 
the sawtooth: h = 6 mm, a — 0.25. The system width is 
10w — 60 mm, one run lasted for 30 internal seconds. The 
plots for n — 0.35, n — 0.7, n = 1.05, n = 1.75, n — 2.8, and 
n = 4.2 are average of 30, 18, 12, 7, 7, 7 runs, respectively. 
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C. Frequency dependence 



Another interesting question is what happens if the 
frequency is varied. A sawtooth shape producing nega- 
tive transport for few particles was chosen. Then with 
different layer widths the horizontal transport shows in- 
teresting behavior as a function of frequency. The results 
are shown in Fig. |Io|. 

Generally, the negative transport becomes stronger 
with increasing frequency. In three cases (n = 1.05, 
n = 1.75, and n — 2.8) the transport reverses, from pos- 
itive to negative. This phenomenon is again very illumi- 
nating when we try to explain the mechanism of negative 
transport. 



D. Dependence on the coefficient of restitution 
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FIG. 12. Horizontal velocity (v x ) as a function of V. The 
shape of a sawtooth: h — 3 mm, a = 0.1. The layer width 
was n — 0.7, the system width 10m = 60 mm. The simula- 
tions ran for 30 internal seconds. A power low function was 
fitted: F* 7 with r* = (T - r c )/T c . The fitted values are 
T c = 1.31 ± 0.03, and 7 = 0.27 ± 0.05. 



We have demonstrated that the transport properties of 
the granular layer in our system depend on the various 
parameters in a non-trivial way. In particular, increasing 
the layer width, through its effect on the enhanced clus- 
tering of the beads, resulted in a change of the direction 
of the current. Thus, we are motivated to investigate the 
effect of changing the coefficient of restitution as well. 
Except for these series, e is fixed at 0.8, but now it varies 
from 0.4 to 1.0. It can be seen in Fig. [ll] that as the dis- 
sipation is decreased the direction of transport reverses 
from positive to negative. The magnitude of the velocity 
reaches a local maximum near e = 0.95, and decreases 
until 1.0. 
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FIG. 11. The horizontal transport as a function of coeffi- 
cient of restitution for three different driving freqencies. The 
layer width is n — 2.8, the shape of the sawtooth and the 
system width are the same as in Fig. [To| The data shown 
here are averaged over 20 runs. 




FIG. 13. Horizontal transport velocity as a funcion of 
driving frequency. The shape of a sawtooth: h = 6 mm, 
a — 0.25. The layer width is n = 1.75. The system width is 
lOiu = 60 mm, the simulations ran for 30 internal seconds. 



E. Phase transition 

When the driving freqency is not high enough (i.e., the 
dimensionless acceleration T < 1), then (v x ) — 0, since 
the whole granular layer has zero velocity compared with 
the oscillating base (solid phase). However, the transport 
may remain zero for T > 1, until a critical acceleration. 
Increasing the frequency, some of the particles, and later 
all of them arc flying (fluidized phase), therefore, trans- 
port may appear. We investigated this phase transition 
(with control parameter (v x )) in case of two different saw- 
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tooth shapes. First we chose a shape which produces 
positive transport for a few particles. Figure [jj] shows 
the r dependence of the velocity. The transport appears 
when r > T c = 1.31. Above this, a power law of the 
form (v x ) ~ (r — r c ) ' 27 can be fitted to the results. 
Choosing a sawtooth for which the transport velocity 



is negative for small particle number, one expects differ- 
ent behavior. Fig. [13| shows how the velocity depends 
on the frequency in such a case. The velocity becomes 
positive growing linearly above the critical frequency, and 
increasing the frequency, it turns back, and the transport 
reverses. 
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FIG. 14. Momentum density field (upper row) and mass density field (lower row) in case of positive horizontal transport 
{{'Vx) = 3.9 ± 0.3 cm/s, h — 3 mm, a = 0.1, / = 25 Hz, n = 0.7). The arrows represent the average momentum vector in a 
box (with width w/6) starting from the center of the box, the radii of the disks are proportional to the average mass. (The 
momentum and mass carried by a particle whose center is in the box are taken into account when calculating the average in 
the box.) The momentum and the mass are averaged for these spatial boxes and phase frames of the oscillation (with length 
2-7T-/6). In the jth frame, the phase <j> = 2nft of oscillation Asm{2n ft) is j2-7r/6 < <f> < { j + l)27r/6. (a) j = 4, (b) j = 0, (c) 
j — 2. The actual height of the simulation cell is three times larger than shown here. 



V. DISCUSSION 

According to our studies of a simplified geometrical 
model the following qualitative argument can be used 
to explain the observed current reversal as a function of 
the particle number: There is an intermediate size and 
asymmetry of the teeth for which a single ball falling 
from a range of near- vertical angles bounces back to the 
left (negative direction) in most of the cases. This ef- 
fect is enhanced by rotation, due to friction between the 
ball and the tooth. However, if there are many particles 
present, this mechanism is destroyed, and on average, the 



direction of the motion of particles will become positive 
(the "natural" direction for this geometry); this corre- 
sponds to the usual ratcheting mechanism characterized 
by larger distances traveled by the particles along the 
smaller slope with occasional jumps over to the next val- 
ley between the teeth p|- |T^ , |23]] . There is no net current 
for symmetric teeth in our case, however the motion of a 
angle particle is very interesting in that situation as well 

The reversal of the current as a function of the fre- 
quency can be interpreted in a similar manner. For 
smaller frequencies there are many particles close to the 
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base and the current is positive. For larger frequencies 
the granular state is highly fluidized, the density con- 
siderably decays and the mechanism leading to nega- 
tive transport for a small number of particles (discussed 
above) comes into play. 

To obtain a deeper insight into the process of horizon- 
tal flow, the momentum density field and mass distribu- 
tion is plotted in case of positive (Fig. [l4| ) and negative 
(Fig. |l5|) transport. In Fig. [l4| it can be seen that the 



particles slide down the side with the smaller slope, and 
when the base is in rising phase, they are given a posi- 
tive net horizontal impulse. It is worth noting that most 
particles do not collide with the opposite side. On the 
other hand, if the sawtooth is higher and more symmet- 
ric, the off-bouncing particles collide with the opposite 
side of the sawtooth, therefore their horizontal impulse 
is reveresed, and the net current becomes negative (see 
Fig. 0). 





(a) 



(b) 



(c) 



FIG. 15. The same as in Fig. [u| but in case of negative horizontal transport {{v x ) = —3.2 ± 0.3 cm/s, h = 6 mm, a = 0.25, 
/ = 30 Hz, n = 0.35). (a) j — 4, (b) j = 0, (c) j = 1. The ratio of the absolute value of the average momentum in a box and 
the length of the arrow is the same as in Fig. |l4[ 



We interpret the existence of maxima in Fig. as a 
function of the particle number as follows: If only a few 
particles are are present, their motion is erratic, with 
large jumps in random directions. As the number of par- 
ticles is increased, due to a inelastic collapse-like process, 
the particles start to move coherently, and this more or- 
dered motion, together with the right frequency, seems 
to give rise to a kind of resonant behavior as far as trans- 
port is concerned. Similar maxima can also be observed 
in molecular motor calculations and simulations. Thick 
layers move more slowly because of inelastic damping. 

The remarkable result that emerges from both experi- 
ment and simulation is that the flow direction can change 



as the layer thickness varies. This is entirely unexpected 
and requires further investigation; the only related be- 
havior of which we are aware is the alternating current 
direction in a model of collectively moving interacting 
Brownian particles in a "flashing" ratchet potential [|l3| . 

In conclusion, we have investigated granular transport 
in a system inspired by models of molecular motors and 
have observed, both experimentally and numerically, that 
the behavior depends in a complex manner on the param- 
eters characterizing the system. These results ought to 
stimulate further research into this fascinating class of 
problems. 
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APPENDIX: COLLISIONS AS EVENTS 

In an event driven algorithm the pace of the simulation 
is determined by the time sequences between collisions. 
In our simulation, which uses Lubachevsky's event driven 
algorithm |2(| with extensions, there are four types of col- 
lisions, i.e., a particle may collide with a (i) particle, (ii) 
resting segment, (iii) oscillating segment, and (iv) oscil- 
lating arc. To find the next collision is the heart of the 
matter. One has to calculate the times of the possible col- 
lisions and find the closest one. Let us denote the time 
of the last event by t^E- Now we show the calculation of 
the collision time in full detail in all cases. 



1. Particle— particle collision 

The motion of particle i until its next collision in a 
gravitational field is determined by its position roi and 
velocity voi at time U. As the two particles (1 and 2) 
are moving along parabolic paths, their position at time 
t are 



(Al) 



ri(i) =r i+voi(f-*i) -jf(<-<i) 2 
r 2 (t) = r 02 +v 02 (<-< 2 ) -j|(<-< 2 ) 2 , (A2) 

where j is a unit vector having opposite direction to that 
of gravity. A collision occurs if Wjjt ) — r 2 (i)| = R\ + i? 2 , 
which can be written, using ( Al - |A2| ) and introducing the 
following notations 

P = voi - v 02 +jg(h - i 2 ) 

q = roi r 02 + v 02 i 2 - vojii - j-(tf - t\) 



as 



pi- 



R\ + R2- 



(A3) 



Taking the square of (A3) and divide it by 2 we get 

^ + p . qt+ S!^±M =0 . (A4) 

The two particles can collide with each other only if 
|p| ^ and 



R 




FIG. 16. Particle colliding with a resting segment. The 
particle arrives from the inside of the simulation space and 
collides with the segment between its ending points si and 

S2- 

is nonnegative, otherwise the particles will not collide. If 
D > 0, then the two times when the condition of touching 
is fulfilled are 



icoll 

c l,2 — 



(A6) 



where tf u < tf u . (If D = 0, then if" = tf u , and this 
means that the two particles only touch each other, so we 
may consider that no collision occurs in this case.) Since 
the smaller solution gives the time of the first touching, 
and only collisions in the future should be regarded, the 
actual time of the collision is given by if" if if > thE, 
otherwise no collision occurs. Summarizing the results, 
two particles collide, and then the collision time is if" , 
if and only if all of the following three conditions are 
fulfilled: |p| ^ 0, D > 0, and t\ M > t LE . 



2. Particle— resting segment collision 

Let us denote the position and the velocity of the parti- 
cle by Tq and Vo respectively at time to- Let n be a unit 
vector perpendicular to the segment, pointing towards 
the inside of the simulation space, and Si and s 2 the left 
and right side end point of the segment with respect to 
n (see Fig. [l6|). Since the particle is moving along a 
parabolic path in the gravitational field, its position as a 
function of time is 



r(t) = r + v (i - t ) 



while its speed is 



v(i) = v - jg(t-t ). 
The particle collides with the line if 
(r(t) - bi) • n = R, 



(A7) 



(A8) 



(A9) 



D= (p-q) 2 -p 2 (q 2 -(i? 1 +i? 2 ) 5 



(A5) 



which, after inserting r(t) from (A7) and using the fol- 
lowing notations: 
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-j • n, v= (jgt + v ) ■ n, 



w = (r - si - v < - j|*o) • n - i?, D = v 2 
can be written as 



ur + «t + w = 0. 



Auv 



(AlO) 



We have to consider several cases. First, if u = (in other 
words j and n are perpendicular to each other, therefore 
the segment is vertical) , then at most one solution is pos- 
sible. If v = in this case, we can say that no collision 
occurs. If v 7^ 0, the solution is 



J.COLL 

''0 — ~ J 



(All) 



and this is valid only if £q > tLE- If it 7^ 0, then, 
provided that D > 0, there are different solutions: 



-v±VD 
2u ' 



(A12) 



If D < 0, then we have either one (the particle just 
touches the line) or no solution. In th ese c ases no col- 
lision occurs. Now let us assume that ( AlO ) has one or 



more solutions, which can be ffj , if , or t?f . Never- 
theless, two more conditions have to be fulfilled at the 
same time: 

(i) at the moment of collision the particle has to come 
from the inside of the simulation space, i.e., 



v(i< 



°")-n<0, 



(A13) 



where t denotes a solution. It is clear that even 



if there are two solutions, after (A13) only one re- 
mains. 

(ii) It is not enough to collide with the line of the seg- 
ment, the particle has to collide with it between its 
ending points si and S2. This is true if 



(r(t 



coll \ 



si)-t>0 and (r(f 



coll \ 



s 2 ) t <0, 
(A14) 



where t = 1 s2 Sl i , a unit vector parallel to the 
s 2 -si|' f 

segment. Summarizing, a particle collides with a 
segment only if (AlO) has solution, and with this 
solution both (i) and (ii) conditions are fulfilled. 



3. Particle—oscillating segment collision 

In this case the time of collision cannot be obtained in 
a closed form, a numerical method has to be used. The 
segment (and so its left ending point) is oscillating in 
vertical direction: 




FIG. 17. A particle can collide with an oscillating segment 
only if it enters the area bounded by the vertical dashed lines 
and its lowest point is between the two horizontal dashed 
lines. 



si(i) = sio + jAsin(27r/t + <j) Q ) 



(A15) 



Let the particle has the same properties as in the pre- 
vious (particle-resting segment collision) case. Then we 
get the equation for collision by replacing S! by Si(t) in 
(|AS|) . We can assume that jn ^ 0, so the segment is not 
vertical (if it is, then with a slight modification it can be 
taken as a resting segment). Using the notations 

9 , (M)+v )-n 
U = -2A> b= 4i-n ' 
(2r - 2s 10 - 2v <o - jfl^o) ' n 

w= ' 

the equation is 

ut 2 + vt + w = A sin(2n ft + <j) ). (A16) 

This equation can be solved only numerically. First, we 
have to find the time region in which a solution is possi- 
ble. In Fig. [n] one can see the segment in the lowest and 
the highest positions, and two particles touching the seg- 
ment at its left end point su (in the lowest postion) and 
right end point s 2u (in the hightest position). It is easy 
to see that the particle can collide with the segment if 

(i) it enters the area bounded by the two vertical dashed 
lines and (ii) the particle's lowest point is in the area 
bounded by the two horizontal dashed lines. Regarding 
these dashed lines as resting segments, using our previous 
results, we can determine the time intervals T v and XK 
when the particle is in a position that conditions (i) and 

(ii) are fulfilled. For this, one point of each of these seg- 
ments has to be determined (as their direction is known) . 
The center of the particle touching the segment at point 
Si; is at point Si;+i?n, so one point of the vertical bound- 
ary line is su + Rn + R(j ■ n)i, and that of the horizontal 
boundary line is su +i?n — R(j ■ n)j. Similarily, for point 
s 2u , these are s 2tl +i?n— i?(j • n)i and s 2u +i?n— i?(j • n)j, 
respectively. The dot product jn is needed to ensure the 
validity of our results if n is pointing downwards. (In 
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this case 'left' and 'right', and 'up' and 'down' also have 
to be swapped.) The time interval T v , when the particle 
is between the vertical dashed lines, is a closed interval 
(in a special case, when the particle is moving vertically, 
this can include all time). The time interval T/j, while 
the lowest point of the particle is between the horizontal 
dashed lines, can be empty or one closed interval or two 
disjunctive closed intervals. Furthermore, only solutions 
in the future are valid, i.e., a solution has to fall into 
the time interval T* = [t^E , oo[. Summarizing, a particle 
collides with a segment if QA16 ) has solution in the time 
interval T = Th H T v PI T*, and if there are more than 
one of them, the smallest one is the time of collision. In 
our simulation the bisection method |m| was used, with 
ending tolerance 10~ 12 s. 



be calculated similarily to the case of particle-oscillating 
segment collision. 




2A 



FIG. 18. A particle can collide with an oscillating circle 
only if it enters the rectangular area bounded by the vertical 
and horizontal dashed lines. 



4. Particle— oscillating circle collision 

The method of finding the collision time is similar to 
that of in the previous case. The center of the circle is 
moving according to 



s(t) = s +j^sin(27r ft + fo) 



(A17) 



(see Fig. |I|). The particle with radius R\ is mo ving again 
along a parabolic path, its position is given by (A.7). The 
equation for the collision time is 



|r(t) — s(t)| = Ri+R 



2- 



(A18) 



where i?2 is the radius of the circle. It is a transcen- 
dental equation, and we use again the bisection method 
for solving it. The collision may occur only if the the 
particle enters the rectangular area bounded by the two 
horizontal and vertical lines in Fig. [l8|. It is easy to de- 
termine these lines: the horizontal lines contain points 
s ; — R2] and s u + R2], while the vertical lines contain 
points So ± R%i. The time region of possible collision can 
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